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Abstract 

We derive an asymptotic formula for the amplitude distribution in a fully nonlin- 
ear shallow-water solitary wave train which is formed as the long-time outcome of the 
initial-value problem for the Su-Gardner (or one-dimensional Green-Naghdi) system. 
Our analysis is based on the properties of the characteristics of the associated Whitham 
modulation system which describes an intermediate "undular bore" stage of the evolu- 
tion. The resulting formula represents a "non-integrable" analogue of the well-known 
semi-classical distribution for the Korteweg-de Vries equation, which is usually ob- 
tained through the inverse scattering transform. Our analytical results are shown to 
agree with the results of direct numerical simulations of the Su-Gardner system. Our 
analysis can be generalised to other weakly dispersive, fully nonlinear systems which 
are not necessarily completely integrable. 



1 Introduction 

It is widely appreciated that, although completely integrable systems may successfully cap- 
ture many important features of nonlinear dispersive wave propagation, they may fail to 
provide an adequate description of large amplitude waves. Consequently, significant efforts 
have recently been directed towards the derivation and analysis of relatively simple models, 
enabling the quantitative description of the propagation of fully nonlinear waves. Although 
such "intermediate" models are typically not integrable using an inverse scattering trans- 
form (1ST), they often provide the possibility of obtaining important particular solutions, as 
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well as having other advantages compared to the full physical system from the viewpoint of 
numerical simulations. 

One such system, the so-called Green-Naghdi (GN) equations, describes large-amplitude 
shallow water waves [1] and also appears in a number of other fluid dynamics contexts, such 
as "continua with memory" solar magnetohydrodynamics [3j, bubbly fluid flows [4] and 
the dynamics of short capillary-gravity waves [5]. The two-layer versions of the GN system 
m E] provide a broad field for the modelling of large amplitude interfacial waves. It should 
be noted that the one- dimensional version of the GN system, which is the main subject of 
this paper, was originally derived by Su and Gardner [9] using a long-wave asymptotic 
expansion of the full Euler equations for irrotational flow (see also El, Grimshaw & Smyth 
[To]), while the original 2D GN system was derived using the "directed fluid sheets" theory, 
which does not require a formal asymptotic expansion, but instead imposes the condition 
that the vertical velocity has only a linear dependence on the vertical {z) coordinate, and 
simultaneously assumes that the horizontal velocity is independent of {z). In this paper 
we shall use the term SG system which appears historically more correct, at least for one- 
dimensional dynamics. The SG system has the form 

Vt + {w)x = , (1) 



1 

Ut + uu^ + r/x = - 
V 



X 



where, in the context of shallow-water waves, rj is the total depth and u is the layer-mean 
horizontal velocity; all variables are non-dimensionalised by their typical values. The first 
equation is the exact equation for conservation of mass and the second equation can be 
regarded as an approximation to the equation for conservation of horizontal momentum. 
The system ([T]) has the typical structure of the well-known Boussinesq-type systems for 
shallow water waves, but differs from them in retaining full nonlinearity in the leading-order 
dispersive term on the right-hand side of ([1]). We stress that there is no limitation on the 
amplitude assumed in the derivation of ([1]). This, along with the appearance of this system in 
different physical contexts, suggests that equations ([1]) represent an important mathematical 
model for understanding general properties of fully nonlinear fluid flows beyond the strict 
shallow- water limit. For this reason, it is instructive to study its solutions for the full range of 
amplitudes, although in the particular context of shallow- water waves the system ([1]), as for 
any layer-mean model, is unable to reproduce the effects of wave overturning and becomes 
nonphysical for amplitudes greater than some critical value. 
The system ([T]) has a solitary wave solution of the form 



r] = r]o + a sech^ ( -^^={x - ct)] , (2) 

where the solitary wave speed c is connected to the amplitude a by the relationship 



c = Uo + y/rjQ + a . (3) 

Here uq and rjo are the background flow parameters. We note that formula ([3]) appears in 
Rayleigh [TT]. 

The asymptotic reduction of the SG system ([1]) for weakly nonlinear waves is obtained 
by using the standard scaling r] = r]o + 5(, u = uq + 5U, x — cot ^ 5~^^'^x, t 5~^^'^t, where 
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(5 is a small parameter and Cq = Uq + -^/Co is the linear longwave speed. For uni-directional 
propagation, this leads to the KdV equation (see, for instance, Johnson [T^ ) 

Ct + ^ccx + kxxx = . (4) 

With the further scaling x' = \/6x,t' = \/6t,v = 3^/2, (j4]) reduces to the standard KdV 
form 

Vt + VV^ + f = , (5) 

on omitting the primes. Then if one considers sufficiently rapidly decaying positive initial 
data for ([5]) 

f (x, 0) = fo(x) > , vq ^ as |x| — s> oo , (6) 

where vo{x) is a continuous "one-hump" function, then as t ^ oo the solution asymptotically 
consists of a finite number of solitons plus dispersive radiation. 

The associated linear spectral problem for ([5]) is given by the Schrodinger equation 

6ij,, + {v{x,t)-X)ij = 0. (7) 

Note that this is in non-standard form due to the factor 6 in the first term. Then the 
parameters of the solitons and the radiation are found from the 1ST using the initial condition 
Vo{x) ([6]) for f in ([7]) . In particular, a soliton with amphtude a„ corresponds to the eigenvalue 
An = an/2 of the Schrodinger operator. 

If one considers a large-scale initial distribution such that e = l/(A^^'^w) <^ 1, where 
A oc maxfo, and w is the typical width of vo{x), then the contribution of the radiation 
is exponentially small in e and as t ^ oo the asymptotic outcome consists only of a large 
(A^ oc e~^) number of solitons, whose amplitudes can be found from the Bohr-Sommerfeld 
semi-classical quantization rule, 

-^^ V^o(x) - A t/x = 7r(n - i) , n= 1,2,3,. ..,iV. (8) 

Here the integral is taken between the turning points xi^2 defined as the roots of the equation 
vo{x) = A. For every given n formula ([8]) yields a bound state An > 0; A Ai > A2 > ■ ■ ■ > 
Atv ^ with the total number of bound states in the potential vq{x). In the KdV context 
the distribution (IHl), through an = 2A„, gives the amplitude of the n-th soliton in the soliton 
train as t — > 00. For the largest amplitude soliton we have the classical result 

amax = 2Ai = 2A . (9) 

On the other hand, for a given A G [0,^4], the formula ([H]) determines the total number of 
bound states n in the spectral interval (A, A) as 

n = i.p. (F(A) + J), F{X) = ^ [ ' ^vo{x)-Xdx. (10) 
^ VDvr Jxi 

Here "i.p. (■)" denotes the integer part. 
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For large n the bound states are located close to each other and we can introduce a 
continuous amplitude function a = 2A. Then differentiating ([H]) with respect to a, one 
obtains the distribution function for soliton amplitudes in the train 



so that f{a)da yields the number of solitons having amplitudes in the interval (a, a+da). The 
total number of solitons in the soliton train is estimated by setting A = in ([8]) (alternatively, 
one can integrate /(a) from to Omax) 



The distribution (ITT]) was originally introduced into soliton theory by Karpman (see 
also Whitham [1A\ and Karpman [I5]). More recently, distributions of this type have been 
obtained, also from the associated linear spectral problem, for the defocusing nonlinear 
Schrodinger equation [12] and for the Kaup-Boussinesq system [T7] . 

One can see that integrability of the KdV equation plays a crucial role in the derivation of 
the formulae fl8l)- f|T2l) through the eigenvalue problem for the associated linear Schrodinger 
equation ([7]). At the same time one can observe that the formulae obtained correspond to 
the semi-classical approximation and that this very approximation can be applied directly to 
the KdV equation, by-passing the 1ST formalism. It is well-known (see Lax and Levermore 
[TS]) that the semi-classical limit of the KdV equation with decaying initial data leads to the 
Whitham modulation system, which can be derived by a direct averaging of the KdV con- 
servation laws over nonlinear wavepackets [19j or from a multiple-scale perturbation analysis 
[20l[2T]. Therefore one could expect that the results (p|)- (fT2|) could be obtained within the 
framework of the modulation theory alone. Indeed, it was shown in Gurevich et al [22] and 
El & Grimshaw [23] that the long-time asymptotics of exact solutions to the KdV- Whitham 
equations with decaying initial data agrees, to first order in t^^, with Karpman's formula 
(|TT|) . Still, modulation solutions in the cited papers rely on the integrability of the KdV 
equation as they employ the presence of Riemann invariants for the associated Whitham 
system [TH], and this latter property is due to the unique finite-gap spectral structure of 
periodic (or, more-generally, quasi-periodic) solutions of the KdV equation 

An important theme in the original work of Whitham [TH] is that, unlike the 1ST, the 
nonlinear modulation approach can be applied to nonintegrable dispersive systems provided 
a minimal structure is present, that is the existence of periodic solutions characterised by 
a certain number of parameters, and the availability of a limited number of conservation 
laws. The resulting modulation system consists of hydro dynamic-type equations and is 
hyperbolic in many cases, which implies that it can be treated using classical methods of 
characteristics theory (see, for instance, Courant & Hilbert [25]). This opens the possibility of 
obtaining some analytical results for non-integrable dispersive wave systems in the framework 
of modulation theory. 

Indeed, it was recently shown in El [2S] that the availability of a number of important 
exact results for the single-phase modulation theory is, in fact, due to certain very general 
properties of modulation systems and is not connected with the presence of the Riemann 
invariant structure. These properties have been used in El, Grimshaw & Smyth [10] to derive 




(11) 




(12) 
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the main parameters of large-amplitude shallow-water undular bores evolving from an initial 
step for the SG system ([T]). In the present paper, we extend the modulation analysis of El 
f2E\ and El, Grimshaw & Smyth ^lOj to obtain an analogue of Karpman's formula fllip and 
its consequences @ and ([T2l) for the SG system for an initial disturbance in depth and 
velocity that decays at infinity. Our approach is based on an assumption, confirmed by 
numerical solutions, that the main qualitative features of the KdV evolution of large-scale 
disturbances, such as the formation of a single-phase undular bore and its further evolution 
into a solitary wave train with a negligibly small contribution of the linear radiation, remain 
present in the SG model. Our analytical results will be supported by comparison with full 
numerical solutions of the SG system. 

2 Asymptotic formula for the KdV equation derived 
from modulation theory 

2.1 Characteristic integrals of the modulation system 

It is instructive to start with a demonstration of our general approach to an asymptotic 
description of solitary wave trains by using the KdV equation as a "test" example. Our 
ultimate goal here is to derive Karpman's formula ( iTTl) without invoking the integrability 
properties of the KdV equation and its modulation system. This derivation will then serve 
as a prototype for calculations for the SG system, which is genuinely non-integrable. 

We consider the initial- value problem ([5]) and ([6]). Without being too restrictive, we also 
assume that Vo{x) = for x > 0, which will simplify the following analysis. As before, we 
consider a large-scale initial disturbance, so that (A^^'^w)^^ <^ 1, where A = max (vq) and w 
is the characteristic width of the initial hump. 

The wave evolution from ([5]) and ([6]) leads to wave breaking at some t = tb, which can 
be taken to be t;, = without loss of generality, which is then resolved by an undular 
bore (which is also often called a dispersive shock wave), which is an expanding nonlinear 
modulated wave train with a distinctive spatial structure. Near the leading edge of the 
undular bore the oscillations appear to be close to successive solitary waves, while in the 
vicinity of the trailing edge they are nearly linear. An asymptotic similarity modulation 
solution for the undular bore evolving from an initial distribution in the form of a sharp 
step was first obtained by Gurevich and Pitaevskii (GP) [27] using the Whitham method 
of averaging over periodic nonlinear wavetrains [TH |T9] and then rigorously recovered in 
the framework of the semiclassical 1ST formalism by Lax, Levermore and Venakides (see 
their review [IH] and references therein). In contrast to the Lax-Levermore- Venakides global 
construction, the direct modulation approach of GP does not rely on the 1ST and thus 
has the advantage of potential applicability to non-integrable systems. This advantage was 
realised in El [26] for the simplest class of problems with step-like initial conditions where 
the modulation solution for the undular bore represents an expansion fan, similar to the GP 
solution. 

For the case of positive initial data, which decay at infinity, the whole initial disturbance 
decomposes, at very large times t ^ e~^, into a chain of solitons with a certain amplitude 
distribution. As was mentioned in the Introduction, the semi-classical 1ST approach leads 
to Karpman's formula flTT]) for this distribution. However, it transpired that the same result 
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can be obtained using the derivation of the intermediate modulation solution for the undular 
bore and then considering its long-time asymptotic behaviour [2^ 1^ . In this connection it 
should be mentioned that the modulation dynamics in these problems with decaying initial 
data is not self-similar, so the original Gurevich-Pitaevskii [27] theory required significant 
development in order to be applied to such problems. It transpires that further development 
is required to extend modulation analysis to non-integrable systems with decaying initial 
data. Below we present the GP formulation for the KdV equation in a form convenient for 
further generalistion to non-integrable systems. 

The local waveform of the undular bore is described by the single-phase periodic solution 
of the KdV equation travelling with constant velocity c, that is v{x,t) = v{6), 9 = x — ct. 
This solution is specified by the ordinary differential equation 

{vef = Q{v) , Q{v) = ^{v - v^){v - V2){V3 - v) , v{e + 2Tx/k) = v{e) , (13) 

V:i > V2 > vi being constants of integration. The phase velocity c, the amplitude a, the 
wavenumber k and the mean v are expressed in terms of the polynomial roots Vj as 

f dv \ _ TT /■ vdv ^, 

while the wave frequency is uj{v, k, a) = kc. When a — the solution of (fT3|) turns into 
small amplitude harmonic waves with the dispersion relation 

uJo{v, k) = uj{v, fc, 0) = kv — k^ . (15) 

In an opposite extreme, when — 0, the travelling wave u{6) transforms into a soliton with 
its velocity depending on the amplitude as 

Csoi = ^ + 3 (16) 

If one allows slow dependence of Vi, V2 and or, equivalently, v, k, and a on x and t, 
one arrives at the modulation system which can be derived by averaging two of the KdV 
conservation laws over the periodic wave family (1131) 



and then closing (fT7|) with the equation for "conservation of waves" (see Whitham [Ti] ) 

dk du(y,k,a) 

dx ^^^^ 

We note that the detailed expressions for the averages in ( |T71) will not be needed. To describe 
the evolution of modulations in an undular bore the system ( |T71) and ( |T8l) must be equipped 
with matching conditions ensuring continuity of the mean u at the undular bore edges x^(t) 

x = x-{t): a = 0, v = f3^{t), , . 

x = x+{t): k = 0, v = (3+{t). 
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Here = P{x'^(t),t), where P{x,t) is the solution of the Hopf equation 



A + /?/5x = 0, P{x,0) = vo{x), (20) 

which represents the dispersionless hmit of the KdV equation (jlj) and is vahd outside the 
oscillatory region. As we have assumed that ^0(3;) = for x > 0, we have f5~^ = 0. 

The trailing x~(t) and the leading x~^{t) edges of the undular bore represent free bound- 
aries defined by the kinematic conditions 

^ = c,(/5-, fc-) , ^ = c.,K/5+, a^) , (21) 

where Cg = dujo/dk is the group velocity of linear waves and k~ = k^{t) = k{x^{t),t), 
a"*" = a'^{t) = a{x^{t),t). It is important to note that, according to general properties 
of quasi-linear hyperbolic systems (see, for instance, Courant & Hilbert [25]), the curves 
x^(t), being the lines which match two analytically different solutions, must coincide with 
characteristics of the modulation system f|T7|l and f|T8|) . Indeed, as we shall see, the kinematic 
conditions fl2Tl) and the requirement that the undular bore boundaries must be characteristics 
of the modulation system are equivalent. 

We now consider the modulation equations in two distinguished limits: as a ^ and 
as k ^ — corresponding to the wave regimes at the trailing and the leading edges of the 
undular bore respectively. When a ^ the oscillations do not contribute to the averaging, 
so F{v) = F{v), where F{v) is an arbitrary function, and the modulation system must 
reduce to (see El 1261 for details) 



a = 0, vt + v ■Vx = 0, kt+{uJo{v,k))x = 0. (22) 

The system (1221) has two families of characteristics, dx/dt = v and dx/dt = duoo/dk. The 
first family is consistent with the characteristics of the "external" Hopf equation (1201) which 
transfers initial data from the x axis into the undular bore region in the (x, t)-plane, while 
one of the characteristics of the second family specifies the trailing edge of the undular 
bore (see the first kinematic condition fl2T|) ). Then, according to the general properties for 
the prescription of Cauchy data on characteristics (see, for instance, Whitham [H]), one 
cannot specify the values of v and k at the trailing edge independently. Of course, a similar 
statement is true for the leading edge as well, and this is why the matching conditions ( fT9l) 
are sufficient to determine the evolution of an undular bore regardless of the fact that they 
involve less variables than the modulation system (fTTll . The admissible combinations of the 
values of v and on a characteristic with a = are found by a substitution oik = k{v) into 
(!22l) . which leads to the ordinary differential equation 

dk_ ^ dujp/dv 
dv V — dujQ I dk ' 

Substituting the linear dispersion relation ( |T5i) into ( l23l) one readily obtains the characteristic 
integral 

k{v) = ^l{v-X,), (24) 

where Ai is an arbitrary constant. 
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We now consider the soliton limit — > 0. Here the wavelength 2TT/k tends to infinity, so 
that the contribution of oscillations to the mean value v vanishes, and, similarly to the case 
of vanishing amplitude, we have F{v) = F{v). Hence, we arrive, again, at the Hopf equation 
for V, 

k = 0: vt + v.v^ = 0. (25) 

We finally pass to the limit as A; ^ in the wave conservation law (fT8|) . This limiting 
transition, unlike that as a ^ 0, is a singular one, so that it requires a more careful treatment. 
Firstly we note that the wave conservation law is satisfied identically for = 0, so we need 
to take into account higher order terms in the expansion of ( ITSl) for small k. It is then 
convenient to introduce a "conjugate wave number" (cf. f|T4|) ) 

k = n{ I (26) 




instead of the amplitude a and the ratio A = k/k instead of the original wave number k, so 
that the parameters {v, A, k) form a new set of modulation variables which is more convenient 
for the consideration of the vicinity of the soliton edge of an undular bore than our original 
set {v, k, a). The variable k can be considered as the wavenumber of a "conjugate travelling 
wave" specified by the equation 

Z] =-Q{v), v{e + 27T/k)=vie), ue[v,,V2], (27) 
d6 / 

where Q{v) is the same as in Eq. (fT3l) and 6 = x — ct is a new phase variable characterised 
by the same phase velocity c. Equation (1271) specifies periodic solutions of the "conjugate" 
KdV equation 

Vi + vvi - Wiii = , (28) 

which can be obtained from ([5]) by the change of variables x ^ ix and t it. It is not 
difficult to infer from fl26|) and fl27|) that in the limit Vi V2 (i.e. k 0) k has the meaning 
of an inverse soliton width (or "soliton wavenumber") defined by the asymptotic behaviour 
in the soliton tails v — v ^ exp(— A;|^|) as \6\ ^ 1. Moreover, it follows from fl27j) that 
the dependence of the soliton speed, which coincides with the value of the phase velocity 
c evaluated in the limit as Vi —>■ V2 on its inverse width (and, therefore, from flTBl) . on the 
amplitude a), follows from the conjugate linear dispersion relation uj = ujo{k). The latter 
is obtained from the linear dispersion relation flTSl) by the change of variables k ^ ik and 
ujq —iojQ, i.e. iuJo{k) = uo{ik). We thus have 

UQ(k) = kv + k^ (29) 

and, comparing with fll6p . we obtain the well known relationship between the KdV soliton 
amplitude and its inverse width 

k = 0: ^ = • (30) 

We are now ready to study^the asymptotic expansion of the wave conservation law for A; <^ 1. 
First we substitute k = Ak into Eq. flTSj) to obtain the equivalent representation 

kA^ + uA^ + A{k, + u^) = , (31) 
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where uj = ck. Next we consider Eq. (131]) for small A <^ 1 and assume that A ^ ^t,-^x 
for the solutions of interest. We note that this is known to be the case for modulation 
solutions describing undular bores in weakly dispersive systems, where at the soliton edge 
one has k —>■ 0, but \kx\, \ kz\ oo (see El [26] for a general discussion of this behaviour and 
Gurevich and Pitaevskii [27| for the detailed calculations in the KdV case). Then to leading 
order we obtain the characteristic equation 

OA UodA 

ar + f&="' '''' 



or 



A = Ao on ^ = ^^, (33) 

where Aq ^ 1 is a constant. In particular, when Aq = (that is = 0), the characteristic 
fl33]) specifies the leading edge of the undular bore (cf. fl2T]) ). Now, considering a restriction 
of equation ( |3T]) to the characteristic family dx/dz = ooo/k and using that uj = ujq to leading 
order, we obtain 

r /— N dx u}n(v, k) 

k = 0: kt + iuo)x = on = . (34) 

at k 

We note that the equation kt + {u)q)x = arises as a "soliton wavenumber" conservation law 
in the traditional perturbation theory for a single soliton (see, for instance, Grimshaw [29]). 
but to be consistent with full modulation theory it should be considered along the soliton 
path dx/dt = Csoi = ooo/k. 

Since v and k cannot be specified independently on a characteristic, there should exist 
a local relationship k(v) consistent with the system (120]) and (IM]) . Substituting k = k{v) 
into (IM]) and using (120]) we obtain an equation for k similar to equation (123]) for k obtained 
earlier in the opposite limit as a ^ 

dk dujo/dv 

— = ^. (35) 

dv V- duo/dk 

Substituting ujq (129]) into (135]) one readily integrates to obtain 



k{v) = ^liX,-v), (36) 

where A2 is an arbitrary constant. 

Now we use the fact that the linear wave packet at the trailing edge and the lead solitary 
wave are not independent, but rather are constrained by the condition of being parts of 
the same undular bore. So, if one considers a pair of integrals and (j5B]) with certain 
constants Ai and A2 in the context of the same modulation solution, then Ai and A2 cannot 
be set independently. One can see from that if = Ai, then k = 0, which then must be 
consistent with equation (136]) derived for the soliton configuration. But, at the same time, 
equation (}21]) corresponds to a = 0, which, together with k = 0, implies k = (see (j3U]) ). So 
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setting f = Ai in (136|) immediately yields Ai = A2 = A and we therefore arrive at the set of 
two consistent characteristic integrals of the modulation system 

/l = ^ = 0, k=\r-{v-X)} on ^ = ^, (37) 




k = 0, k=^l=.{\-v)} on Zl = ^^^^l^. (38) 

One should add that both k and k are required to be real, so the integrals /i^2 involve different 
parts of the domain of the function v(x, t). 

2.2 Total number of solitons 

We start with the determination of the total number of solitons generated by the decay of 
the given initial disturbance ([5]) and then proceed by obtaining a more detailed description 
of a soliton train by finding the distribution function for the soliton amplitudes in terms of 
the initial data. 

We consider the wave conservation law (|T8l) . For the case of the decaying initial profile ([5]) 
one obviously has a; — > as |a;| ^ cxd for all t, and hence equation (|T8l) implies conservation 
of the total number of oscillations (wave crests) 



hoo 



N ^ — I kdx = constant . (39) 



We use an approximate equality sign here due to the inherent asymptotic character of mod- 
ulation theory. Next we note that, qualitatively, the process of soliton generation during 
the evolution of the large-scale, decaying initial profile can be described as follows (see 
Gurevich, Krylov & Mazur [30] for a quantitative justification): each wave crest is generated 
at the trailing edge as a vanishingly small amplitude linear wave and asymptotically, as 
t 00, transforms into a soliton. Formula ( |39l) then can be used for the evaluation of the 
total number of solitons in the eventual soliton train. 

To evaluate the integral in fl39|) one needs to know the function k{x,t) for all x at any 
particular t, say at t = 0. The difficulty with the traditional Gurevich-Pitaevskii approach to 
the undular bore description is that the wave number k in this approach is defined only within 
the undular bore region [a;~,a;+]. However, one can extend the notion of the wavenumber 
to the entire x-axis by defining the function k{x,t) in such a way that its behaviour outside 
the undular bore region is consistent with the prescribed values of k along the leading and 
trailing edges x^{t) for all t. 

At the leading edge x^{t) we have k = for all t. Therefore we define k = for x > x^{t). 
Hence here we have 

k{x, 0) = for x> 0. (40) 

At the trailing edge x~{t) the amplitude of the wave a = 0, and since the trailing edge 
represents a characteristic of the modulation system, the value of the wavenumber k~ = 
k{x~ ,t) is determined by the boundary value of v{x~ ,t) = f3~ from the characteristic integral 
h in (1371) for a certain A (see ( ITQil and (120!) for the definitions of P'^it)). To be consistent 
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with the matching conditions (fT9|) . we need to set A = 0. Indeed since at the leading edge 
v{x~^,t) = /3+ = 0, then the dependence k~{j3~) must correctly reproduce the value = 
for the case when (3~ = (5^ = 0, implying a = and k = simultaneously, which by fl3Up 

also implies A; = 0. So we obtain 

k- = ^WW- (41) 

Now upstream of the undular bore, that is x < x~, the function t) satisfies the Hopf 
equation (!20|) with the initial condition 0) = fo(a;), i.e. implicitly P{x,t) = vo{x — [3t). 
Therefore for x < x~{t) we need to define the wavenumber as k{x,t) = k~{f3{x,t)), where 
(3{x,t) is the aforementioned simple-wave solution. This extension basically implies that 
we assume that upstream of the undular bore, where a = 0, the relationship (P7|) holds 
everywhere (not only for a special family of characteristics dx/dt = duj^dk). Indeed it can be 
readily seen that k = vWs) f is a solution of the reduced modulation system (l22l) regardless 
of the particular characteristic family. Thus at t = we obtain 



k{x,0) = k-{f3{x,0)) = \l^vo{x) for x<0. (42) 



We note that since fo(a; > 0) = this formula is also consistent with our definition ( HOl) 
for k{x, 0) for X > 0, so that P2|) gives the function k{x, 0) = kQ{x) = ^/2/3{vo{x)) on the 
entire real line. Now substituting (142!) into (139|1 we obtain 

1 p+co -I r+oo 

N=— ko{x)dx = ^^ vl^\x)dx, (43) 
2vr J-oo V67r J^oo 

which agrees with the 1ST result f|T2|) . 

2.3 Soliton-amplitude distribution function 

We note that the formula fH51) in fact represents a particular case of a more general expression, 
which is obtained by retaining the parameter A in the modulation integral Ji (1371) so that 
instead of ko^x) in (H3l) one introduces 



fco(x. A) = a/ {2/3){vo{x) - A) for xG[xi,X2], 
and A;o(x, A) = for x ^ [xi,X2], (44) 

where X\ < X2 are the roots of the equation fo(x) = A. Thus we arrive at the continuous 
family of quantities characterised by the parameter A G [0, A] 



1 



X2 



F{X) ^ — / ko{x, \)dx = / ^/vo{x) - X dx , (45) 
271" J-oo V67r 

To clarify the meaning of the quantity F{X) in the framework of modulation theory we 
introduce the notion of a A-section of the initial profile, which is simply a segment of the 
function Vq{x) for which fo(x) > A (see Fig. [T]), and consider the integral 



F(A,t) = -^ / V(2/3)(t;(x,t)-A)dx (46) 
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X 




Figure 1: A-section of the initial profile (below) and the behaviour of the corresponding 
confining characteristics for part of the KdV modulation solution (above). r(A) is the domain 
of influence of the interval [xi, X2] 

for t < t2, where ^2 is the time at which the characteristic x = X2 + Xt intersects the trailing 
edge x^(t) of the undular bore x^(t2) = X2 + \t2 (see Fig. [T]). Note that this ^2 plays for the 
A-section of the (evolving) profile v{x, t) a role similar to the breaking time, as the nonlinear 
oscillations emerge in the vicinity of the point {x~{t2),t2) and the system fl22l) ceases to be 
valid. However, unlike at the breaking point, there is no gradient catastrophe at (a;~(t2), ^2)- 
We now see that for t < t2 one has dF/dt = and thus F(A, t) = F(A, 0) = F{X). Indeed 
the integrand in fH^ is nothing but the exact solution of the modulation system fl2^ . 
But this system conserves the integral k(y{x,t))dx, so that we have 

— / k{v{x,t))dx = — / k{v{x,t))dx = — / k{vo{x))dx = F{X) . (47) 
271" J-00 27r J^^j^xt 27r J^^ 

It follows from (1471) that F{X) is the number of linear modulated waves (i.e. the number of 
wave crests) "contained" in the A-section of the initial profile. All these waves are "released" 
into the undular bore during the time interval (^25^1) (see Fig. [1]) and eventually transform 
into solitons as t ^ 00. 

Now we establish a correspondence between the A-section of the initial profile and a 
certain part of the solitary wavetrain as t — > 00. This correspondence follows from the 
detailed consideration of the behaviour of the characteristics for the modulation solution 
made in Gurevich, Krylov & Mazur [30] where it was shown that all the characteristics of 
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the modulation system issuing from the points of the interval [xi,X2] are confined to the 
area r(A) of (x, t)-plane enclosed by the leading edge x~^(t) and the characteristic emerging 
from Xi (see Fig. [T]). On the other hand, it was shown that for sufficiently large t the 
characteristics emerging from the points x ^ [xi,X2] lie entirely outside r(A). Therefore all 
wave crests generated at the trailing edge between the points = X2 + Xt2 and x^ = Xi + Xti 
remain within the described region r(A) for all t and eventually transform into the portion of 
solitons having their amplitudes in the interval [ci(A); a(y4)]. In other words, r(A) represents 
the domain of influence of the interval [a;i,a;2] and thus, formula (H5l) defines the number of 
solitons in this region as t — > oo. Indeed when A = it coincides with formula (l43ll for the 
total number of solitons. 

We still need to find the dependence a(A) to obtain the distribution of solitons as a 
function of amplitude. In this regard we take advantage of the fact that due to the condition 
v{x, 0) = for a; > 0, all solitons propagate on a zero background, so substituting v = into 
I2 we obtain A = 3/c^/2. Now using the connection (130|) between the inverse width of a KdV 
soliton and its amplitude we obtain the relationship a = 2A. In particular we immediately 
recover the 1ST formula amax = cl{A) = 2A (see (^). Then differentiating F{\) with respect 
to a we obtain Karpman's formula ffTTl) for the soliton amplitude distribution function 



/(a) 



dF{a/2) 



da 



1 1 dx 
/ — 

47ry6L ^/vo{x) - a/ 2 



Summarising, we have managed to reproduce the semi-classical 1ST results pertaining to 
soliton dynamics using, technically, only the linear dispersion relation and the characteristic 
velocity of the KdV equation in the dispersionless limit. Of course, we took advantage of our 
knowledge of the qualitative behaviour of the characteristics of the modulation equations for 
the solutions considered here. Next assuming the same qualitative behaviour of modulations 
as a plausible assumption, we extend the above reasoning to the non-integrable SG model ([1]) 
with the aim of obtaining the counterparts of the formulae and fHHl) for fully nonlinear 
shallow water dynamics. Our analjd^ic results will then be compared with full numerical 
solutions of the SG system. 



3 Asymptotic description of a solitary wavetrain in the 
Su-Gardner system 

3.1 Conserved quantities in the SG solitary wavetrain 

We consider the SG system ^ with the initial depth profile 

77(x,0) = rjQ^x) > 1 : 770(2;) ^1 as |x| — > 00 (49) 
and the velocity profile connected with (150|) by the simple wave relationship 

u{x,0)=Uo{x) = 2{^/^-l), (50) 

so that ^0(2;) ^ as |x| 00. Let max[r7o(a;)] = 1 + A and let w be the typical width of 
i]o{x). Similarly to the KdV case we shall assume that A^^'^w ^ 1. Also, for convenience. 
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Figure 2: Solitary wavetrain developing in the SG system from the initial disturbance ([55 
with A = 1 and w = 10. The plot is ri{x, t) at t = 400. 



we assume that rjoix) = 1 for x > 0, which will guarantee that solitary waves propagate on 
the undisturbed background rj = 1. 

The system ([1]) conserves the four physical quantities 



mass 



J, 



u + -rfu^^]dx] 





J. 



[rju]dx] 



irrotationality 



momentum 



■h= I \-[JI — ij" + -^M" + -V"u"^\<ix energy. 
We now consider the SG system with the specific initial conditions 

2(v/^-l 



r7o(a;) = 1 + ^4 sech [x/ty], 



Uq^X) 



(51) 



(52) 



(53) 



(54) 



(55) 



with the aim of verifying numerically some of our main assumptions about the qualitative 
features of the asymptotic behaviour as t ^ oo. A typical long-time outcome is shown in 
Fig. [21 Comparison of the values of the four conserved quantities (I5T1)-( 15^ computed for 
two initial profiles having the form (l55l) with A = 0.4, w = 13 and A = 2, w = 10 with 
their values calculated for the solitary wavetrain at large t shows that the relative change, 
due to radiation, is 0(10^^), which is within the accuracy of modulation theory, valid for 



(A 



1/2 



W 



0.12 and 0.07 respectively. Therefore we can assume that with good 



accuracy a positive large-scale initial disturbance is completely transformed as t — >■ cxo into 
a solitary wavetrain and we can neglect the contribution of the radiation component in the 
modulation analysis. 
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3.2 Modulation characteristic integrals 

We now apply the procedure for the derivation of the amphtude distribution function, de- 
veloped in Section 2 in the context of the KdV equation, to the SG system. The key objects 
are the modulation characteristic integrals, having in the KdV case the form fl37|) and fl38l) . 
The ordinary differential equations for these integrals for the SG case (the counterparts of 
the KdV case equations fl2^ and (I35D) were derived in El, Grimshaw & Smyth [TU], so here 
we present only the resulting formulae. It is important to mention that the derivation in 
El, Grimshaw & Smyth [10] was performed in the context of simple undular bores, so that 
it incorporated the simple wave relationships between the depth and velocity jumps across 
the bore. In the present study this condition is automatically satisfied by choosing initial 
conditions in the form (l49l) and (!50l) . 

The key ingredient of the modulation characteristic integrals is the linear dispersion 
relation for modulations, which in the SG case has the form 

/ ^1/2 \ 

a = 0: u = ujoiv, u, k) = k + ^ ^^^^^^ j . (56) 

Here, as for v in ( ITSl) . r/ and u denote the values averaged over the period of the travelling 
wave. Next we incorporate the simple wave relation u{r]) = 2(a/^— 1) into (1561) to obtain 
the dispersion relation for a linear dispersive wave propagating on a slowly varying simple 
wave background 

Qoiv, k) = MvMv), k) = 2k{t'^ - 1) + ^^^^^l^/^y/^ ■ (57) 

Next, using reasoning identical to that described in Section 2 (see [101 [261 [28] for additional 
details pertaining to bidirectional systems), we obtain the characteristic integrals of the 
modulation system defined by the ordinary differential equations 

^ dk dilo/df] dx d^lo 

- vifj) - dno/dk dt"^' ^^^^ 

fc = 0, f= ^^''^ ^ on ^ = a. (59) 
dv V{f])-dno/dk dt k 

Here 

V{r]) = u{r]) + 77^/2 = 3^1/2 _ 2 (60) 

is the characteristic velocity of the right-propagating simple wave of the ideal shallow water 
equations (i.e. the dispersionless limit of the SG system) and 

~ ~ _ _ '^771/2 

no(f], k) = -inoif], ik) = 2k(f/^ - 1) L (61) 

is the SG "solitary wave dispersion relation", an analogue of the KdV formula (129!) . k being 
the soliton wavenumber. 
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Substituting (1571) into (158|) and introducing a = (l + A;^?7^/3) as a new variable instead 
of k we obtain an ordinary differential equation in separated form 

d5 _ + " + . ,62) 

7] a{l + «)(« — 4) 

Integrating (162|) we obtain 

- 1 //I \ 21/10 , \ 2/5 



Ai 



(63) 



wfiere Ai is an arbitrary constant of integration. Similarly, substituting (l6Tl) into (l59l) and 
tlien using 5 = (1 — /c^rJ^/3)~^/^ instead of A; we arrive at the same separated ordinary 
differential equation (!62|) . but now for for 5(7j) 



*; 2(1 + 3 + 3^) i^^j 



with the integral 



rj a{l + a){a — 4) 

_ 1 //I ~\ 21/10 , ~N 2/5 

A2 5^ 



V2I 3 J I 2 J 



(65) 



A2 being another constant of integration. 

Now we recall that the constants Ai and A2 cannot be set independently if the character- 
istic integrals (!63l) and (!65|) are considered in the context of the same modulation solution. 
One can see from the zero amplitude integral I2 (163!) that r/ = Ai implies a = 1 and, there- 
fore, k = 0. But = 0, in its turn, corresponds to the solitary wave limit, so the equality 
rJ = Ai must be consistent with the zero amplitude reduction of the "soliton" integral I2 
(|65|) . To obtain this reduction we calculate the solitary wave velocity using the conjugate 
dispersion relation fl6T]) 

c=^M2^ = 2(ryV2_i)_^i/2~_ (gg) 
k 

On the other hand, we have from (I3l), after replacing rjo with r] and uq with u = 2{rf^'^ — 1), 



c = 2{rf'^ + ^r] + a. (67) 

Comparing (l66l) and ( l67j) we obtain the expression for the solitary wave amplitude in terms 
of the variable 5, an analogue of the KdV formula ( |30ll . 



= : a = 77(52 - 1) . (68) 

Now one can see that the reduction as = 0, a = implies 5=1 which, by fl65l) . immediately 
yields = A2. Therefore, similar to the KdV case, we have Ai = A2 = A and the relationships 
(|63|) and fl65l) become a set of two consistent modulation characteristic integrals parametrised 
by the same value A (cf. fl37|) and fl38|) for the KdV case). 
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3.3 Total number of solitary waves 

As we described in Section 2.2 the total number of solitary waves in the soliton train evolving 
out of the initial large-scale disturbance fj^9l) and fl50|l can be found from the modulation 
formula 

^ p+oo 

N^— ko{x)dx, (69) 
2vr J-oo 

where ko{x) is the initial distribution of the wavenumber corresponding to the initial condi- 
tions ( l49i) and ( l50i) for the depth rj and velocity u. This distribution is found from the group 
velocity characteristic integral of the zero amplitude reduction of the modulation equations. 
In the case of the SG system this is the integral Ji given by formula (163|) . in which one 
assumes Ai = A = 1. The latter equality follows from the requirement that the integral Ji 
must correctly reproduce the mean value r/ = 1 for the solitary wave background when one 
sets k = (i.e. a = 1) in (!63ll . Next one replaces fj in Ji with its distribution rio{x) at t = 
to obtain an implicit expression for ko{x) in terms of rio{x) 

T / A \ 21/10 /, , N 2/5 

,„W = ^(i^) 0^) , (70) 



where ao{x) is connected with ko{x) via the relationship 

ao = {l + {kom{x)f/3)-'/\ (71) 
Thus expressing ko from ( 17T1) and substituting it into ( 169|) we obtain 



1 Z"^" \/3(l-a^(x)) , 
27r y„oo r]Q{x)ao{x) 

which, together with (170]) . completely defines the total number of solitary waves forming in 
the initial value problem ([1]), fj^9l) and fl50l) . 

The first term of the expansion of fl72|) and flTOj) in 770(2;) — 1^1 (i.e. 1 — ^ 1 by 
[DD) yields 

3"l 



^-\/2^j Vvo{x)-ldx, (73) 

corresponding to the KdV result fH51) . taking into account that the weakly nonlinear reduc- 
tion of the SG system yields the KdV equation in the form ([5]), so to obtain (HSl) one needs 
firstly to set rjoi^x) — 1 = (o{x) in (1731) and then apply the rescaling as described in Section 1. 

To compare our analytical results with full numerical solutions of the SG system we 
consider the dependence of on A^^'^w, which is the parameter suggested by dimensional 
analysis (see also (fT2l) ). for initial data of the form (l55l) . The results are shown in Fig. 3. 
The comparison shown in the left panel is for the solitary wavetrain developing from the 
initial profile fl55|) . with the width set to if = 10 and the amplitude A varying so that A^^'^w 
varies in the interval [6, 13], while in the right panel we fix the amplitude to A = 0.4 and 
vary the width w so that the quantity A^^'^w has the same range [6, 13]. One can see that 
in both cases the modulation formula predicts the total number of solitary waves very well. 
The over-prediction by just one in the second comparison is clearly within the accuracy of 
the modulation approach, which is 0(^4^/^^)"^). 
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Figure 3: Total number of solitary waves forming due to the decay of the initial disturbance 
(l55l) . The initial parameters are w = 10 and A varying (left); A = 0.4 and w varying (right). 
Solid line is the modulation formula (1721) : symbols are the numerical solution; the dashed 
line is the formula fl73|) corresponding to the KdV approximation (jl]). 



3.4 Amplitude distribution function 

Next, following the procedure of Section 2.3, we can "upgrade" formula fl72|) by using in fl69|) 
the full expression k^^x, A), which is obtained from the modulation characteristic integral Ii 
(15^ by retaining the parameter A and replacing the mean depth rj with its initial distribution 
rio{x). The expected result is the SG analogue of the integrated KdV soliton amplitude 
distribution function (H5|) . 

We thus assume that in (17T!) ao = ao{x, A) and then obtain from (163!) an implicit expres- 
sion 



/ \ 1 / A \ 21/10 /T , \ 2/5 

rjoixj 1 /4-ao\ j I + ao\ 



A ^/ao V 3 / \ 2 / 
Then the integrated distribution function has the form 



(74) 



F(A) ^ — / ko{x,X)dx, (75) 



27r 

where the function ko{x, A) is expressed in terms of ao{x, A) as 

(see( l7TI) ) and xi,2(A) are the roots of the equation ko{x, A) = 0. It immediately follows from 
( 1761) and ( 1741) that these roots coincide with the roots of rio{x) = A. So the range of A is 
[1, 1 + A], where A = max[r7o(a;)] — 1. 

Next we establish the connection between the parameter A and the lower bound of the 
amplitude range in the portion of solitary waves corresponding to the "A-section" of the 
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Figure 4: Comparison of the lead solitary wave amplitude function of the initial 

amplitude A. The initial conditions are fl55l) with w = 15. Solid line is the analytical 
solution ( !78l) . the dashed line is the numerical solution. 

initial profile t]q{x). This is found from the characteristic integral I2 (l65l) . in which we set 
A2 = A and r/ = 1 as here the solitary waves are on a unit background (recall that riQ{x) = 1 
for X > 0). Then substituting 5 = a/1 + a from (l68ll we obtain the required relationship 

/ o \ 21/10 / ^ N 2/5 
X = il + ay/'( ' ] ( , ' ] . (77) 



A-y/l + aJ \l + ^/l + a 

Since max X = 1 + A, the amplitude of the largest solitary wave is found from the 
equation 

/ o \ 21/10 / r, X 2/5 

l + A={l + aJ'/U- ( , , ' , . (78) 



Expanding ( 1781) in Om we obtain to leading order am = 2A, as expected from weakly nonlinear 
KdV theory. The relation ( ITHj) is shown in Fig. H] along with the amplitude curve obtained 
from direct numerical solutions of the SG system ([1]) with the initial conditions (155|) with 
w = 15. One can see that the agreement is very good for solitary wave amplitudes up to about 
a = 1.6. The departure of the analytical curve from the numerical one for larger amplitudes 
will be discussed later. One can also observe that the KdV lead soliton amplitude curve 
am = 2y4 gives a very good approximation to the numerical SG curve, in fact even better 
than the modulation SG formula flTHl) . One should, however, bear in mind that formula fl78|) 
for the lead solitary wave amplitude should be used in conjunction with expressions ([2]) and 
(IH]) defining the SG solitary wave profile and velocity and, as such, then provides a consistent 
approximation of the SG solution (and therefore of the full Euler equations solution). At 
the same time the KdV formula am = 2^4 should be considered together with the familiar 
KdV soliton profile and the speed-amplitude dependence (fT6|l which are known not to give 
very good approximations to the SG solitary wave for large amplitudes. 
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Figure 5: Number A/'(a) of solitary waves with amplitudes in the interval [a, a^] generated 
from the initial conditions ( |55l) with A = 0.4, w = 13 (left) and A = 2, w = 10 (right). Solid 
line is the modulation formula ( |75l) : symbols are the numerical solution; dashed line is the 
integrated Karpman formula fl79l) . 

Next putting ( I75i) -( I77|) together we obtain the formula for the number of solitary waves 
with amplitudes in the interval [ 7V(a) = F(A(a)). Indeed 7V(0) = F(l) = A^, where 

is the total number (!72|) of solitary waves in the train. Expanding (!75l) - (!771) in a for a ^ 1 
we obtain to leading order 



where xi^2{a) are the roots of the equation rio{x) = 1 + a/2. Equation fl79l) is the integrated 
Karpman formula for the KdV equation in the form (jl]) (cf. fHSj) ). Comparisons of the inte- 
grated amplitude distribution function A/'(a) with the numerically found number of solitons 
are shown in Fig. Ofor two different initial profiles (one with A = 0.4, w = 13 (left panel) 
and another one with A = 2, w = 10 (right panel)). One can see that the agreement in 
both cases is quite good (again taking into account the accuracy of the modulation approach 
itself). Moreover, one can see from the comparison in the right panel that the modulation 
formula appears to work reasonably well far beyond the range of formal applicability of the 
GP formulation of the problem used here. Indeed, it was shown in El, Grimshaw and Smyth 
(lO j that starting from some critical depth jump Acr ~ 1-4 across the SG undular bore, 
the oscillatory structure of the bore qualitatively changes due to linear degeneration of the 
characteristic field and formation of a rapidly varying, finite amplitude rear wavefront, as 
opposed to the usual vanishing amplitude trailing wave packet assumed in the GP formula- 
tion. As an estimate we can assume yl/2 = A — 1, the value of an "effective jump" taken at 
the level of half the maximum value of the depth disturbance. Then for solitary waves with 
amplitudes greater than Ucr ~ 0.8 one can expect some discrepancy between the modulation 
predictions and the results of direct numerical simulations, as indeed seen in Fig. HI 




(79) 
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Figure 6: Amplitude profile of the SG solitary wavetrain at t = 1500. The parameters of 
the initial profile are A = 2 and w = 10. Solid line is the formula fl8ip : symbols are the 
numerical solution; dotted line is the KdV soliton train amplitude profile fl83|) 

The solitary wave density distribution function /(a) is obtained by differentiating (jTHD 



so that f{a)da gives the number of solitons with amplitudes in the interval [a, a + da]. If 
we let the number of solitary waves per unit length in the solitary wavetrain be k, then it 
follows from the balance relationship ndx = f{a)da that k = f{a)da/dx. Since the speed of 
an SG solitary wave propagating against a background r] = 1 and m = is c = vT+~a, the 
amplitude profile in the solitary wavetrain for x, t ^ 1 is found from the general formula 



A comparison of formula (IHT!) for the solitary wavetrain amplitude profile with the numerical 
profile for A = 2 and w = 10 ai t = 1500 is shown in Fig. [6] 
Next, for the number of solitons per unit length we obtain 



VTTa = x/t + 0{t~^) as 



a ^ {x/t f - 1 . 



(81) 




fiix/tf-l). 



(82) 
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Figure 7: Spatial density of solitary waves k{x, t) in the solitary wavetrain formed out of 
the initial profile (ESI) with A = 0.4, w = 13 (left, t = 1400) and A = 2, w = 10 (right, 
t = 1500). Solid line is the SG formula (1821) : symbols are the numerical solution; dotted line 
is the KdV formula (1831). 

We note that the weakly nonlinear counterparts of formulae (18T1) and ( l82l) corresponding to 
the KdV equation in the form (jlj) are 

a ^ 2{x/t - 1) , ^/w(2(x/t - 1)) , (83) 

where fKdv{o) is given by the Karpman formula (|TT|) . Comparisons of the curve (l82l) with 
the numerically found spatial density of solitary waves in the SG solitary wavetrain for two 
different initial profiles are shown in Fig. [71 As in the previous comparisons, the agreement 
is excellent for moderate {A = 0.4) and reasonably good for very large {A = 2) amplitudes 
of the initial disturbance ( l55l) . 

4 Conclusions 

We have developed a general method for obtaining an asymptotic description of solitary 
wavetrains for initial value problems for weakly dispersive nonlinear systems which may not 
be integrable via the 1ST. The method is based on the properties of the characteristics of the 
associated modulation (Whitham) system describing an intermediate undular bore stage of 
the evolution. We have demonstrated the effectiveness of the developed approach by firstly 
recovering the semi-classical 1ST results for the soliton amplitude distribution function for 
the KdV equation and then by applying it to the Su-Gardner (SG) system ([1]), the ID ver- 
sion of the Green-Naghdi equations, describing fully nonlinear shallow water waves. The 
SG system is not integrable by the 1ST, but has enough structure (a periodic travelling 
wave solution family and four conservation laws) to be amenable to a nonlinear modulation 
analysis. It also transpires that the SG system represents an important mathematical model 
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for understanding general properties of fully nonlinear fluid flows beyond the strict shallow 
water limit. For this reason we have studied its solutions for the full range of amplitudes, 
although in the particular context of shallow water waves the system ([T]) is unable to repro- 
duce the effects of wave overturning and becomes nonphysical for amplitudes greater than 
some critical value. 

The resulting asymptotic formulae for the solitary wave amplitude distribution function, 
the total number of solitary waves and the amplitude of the leading solitary wave in the 
solitary wavetrain formed as the long time outcome of the decay of an initial large scale, 
one hump depth disturbance have been compared with the results of direct numerical sim- 
ulations of the Su-Gardner system. Very good agreement between the modulation solution 
and numerical results have been demonstrated, even for the range of initial amplitudes sig- 
nificantly exceeding that for the applicability of modulation theory due to formation of a 
rapidly varying wavefront at the trailing edge of the undular bore at an intermediate stage 
of its evolution (see El, Grimshaw and Smyth |10j). 

The results obtained for fully nonlinear shallow water theory have also been compared 
with their weakly nonlinear counterparts which are well known from the semi-classical 1ST 
approach to the KdV equation, where they were obtained as consequences of the famous 
Bohr-Sommerfeld quantisation rule (see Whitham [14j and Karpman |15j). Overall, one can 
conclude that KdV theory gives a very good prediction for the lead solitary wave amplitude 
and the spatial density of solitary waves in the fully nonlinear solitary wavetrain, but con- 
sistently over-predicts the number of solitary waves in a given amplitude interval. Indeed 
this discrepancy grows with the initial amplitude. One can also note that this comparison is 
consistent with the modulation results for SG undular bores in the step resolution problem 
studied in El, Grimshaw and Smyth [10], where it was shown that the SG undular bore 
is noticeably narrower and contains less wavecrests than its KdV counterpart for moderate 
to large values of the depth jump across the bore. Also, due to the essentially different 
amplitude-speed relationships for SG and KdV solitary waves (cf. ([3]) and f|T6l) . the former 
being actually equivalent to that for the full Euler theory [llj), the spatial amplitude profile 
of the SG solitary wavetrain is parabolic, in contrast to the KdV classical triangle profile 
(see Whitham [T3] and Fig. 6). 
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